############ Izhikevich model #############
# input 
#	var: c(v,u)
#	p: params c(a,b,c,d,I)
I.ode.system <- function(t, var, p) {
	v <- var[1]; u <- var[2]; 
	a <- p[1]; b <- p[2]; c <- p[3]; d <- p[4]; I <- p[5]

	# reset condition
	if (v >= 30) {
		v <- c
		u <- u + d
	}

	return(list(c(
		0.04 * v^2 + 5*v + 140 - u + I,
		a * (b*v - u))))
}

I.ode.root <- function(t, var, p) {
	v <- var[1]
	return(v - 30)
}

I.ode.event <- function(t, var, p) {
	v <- var[1]; u <- var[2]
	c <- p[3]; 	d <- p[4]
	v <- c
	u <- u + d
	return(c(v,u))
}

I.nullcline <- function(x, p) {
	v <- x
	b <- p[2]; I <- p[5]

	return(c(
		0.04 * v^2 + 5*v + 140 + I,
		b*v ))
}

########### run experiment ##########

	test.ode.system	= I.ode.system 
	test.ode.root 	= I.ode.root
	test.ode.event	= I.ode.event
	test.nullcline	= I.nullcline
	test.params		= c(0.1,0.25,-65,2,0)
	test.xrange		= c(-80,40)
	test.yrange		= c(-30,30) 
	test.xstep		= 5
	test.ystep		= 2.5
	test.xlabel		= "v"
	test.ylabel		= "u"
	test.x.ini		= c(-65,-55,-45)
	test.y.ini		= c(-17,-17,-17)
	test.times		= seq(1,100,by=0.1)
#	test.times		= seq(1,300,by=0.01)  # for finer detail


draw.phase.plane(
			ode.sys = test.ode.system, ode.event = test.ode.event, ode.root = test.ode.root,
			nullcline = test.nullcline, 
			params = test.params,
			xrange = test.xrange, yrange = test.yrange, 
			xstep = test.xstep, ystep = test.ystep,
			xlabel = test.xlabel, ylabel = test.ylabel,
			x.ini = test.x.ini, y.ini = test.y.ini, times = test.times)

#draw.trajectory.over.time(
#						ode.sys = test.ode.system, ode.event = test.ode.event, ode.root = test.ode.root, 
#						params = test.params,
#						x.ini = test.x.ini, y.ini = test.y.ini, 
#						times = test.times,
#						xlabel = test.xlabel, ylabel = test.ylabel)

